Purpose of this analysis

This report fits bivariate Bayesian model similar to those described in Riley et al (2008), Reitsma et al (2005) and Owen et al (2018) in order to get a meta-analysis estimate of the operational characteristics of diagnostic tests for liver fibrosis and cirrhosis based on different biomarkers.

Below we define 4 models, corresponding to those from Riley et al (2008). Each of those models would be fit several times: for each combination of reference diagnostic and experimental diagnostic.

Reading the data extracted during the systematic review

datFile<-"../2021-06-24_RevisedData/HEPSANET_biomarker_cleaned_data.csv"
dat<-read.csv(datFile)

dat<-dat %>%
  mutate(ReasontotestHBVSCLF=gsub(pattern=" ",replacement="",case_when(ReasontotestHBVSCLF==""~"N",ReasontotestHBVSCLF=="I"~"L",TRUE~ReasontotestHBVSCLF))) %>%
  mutate(
    testReason=factor(ReasontotestHBVSCLF,levels=c("C","S","L","F","N")),
    studyType=factor(hosp_com)
    )

setOfDummyVars<-model.matrix(~studyType+testReason,data=dat)
dat<-cbind(dat,setOfDummyVars[,-1])

The file ../2021-06-24_RevisedData/HEPSANET_biomarker_cleaned_data.csv contains 2990 rows / observations and 91 columns / variables.

Notes (as per Alex’s email from 24 June 2021):

dat$testReason[dat$Nomerged %in% c("H014","H009")]<-"S"

dat<-dat %>%
  mutate(testReason=factor(as.character(testReason),levels=c("C","S","L","N")))

Filtering out participants

We filter out participants who are positive for hepatitis C and/or D infection using the variables antihcv, hcvrna, and antihdv (all participants positive for any of these are filtered out). Participants with no records for these variables are retained.

nBefore<-nrow(dat)

dat<-dat %>%
  filter((is.na(antihcv) | antihcv!=1) &
           (is.na(hcvrna) | hcvrna!=1) &
           (is.na(antihdv) |antihdv!=1))

This has filtered out 45 participants, leaving 2945 participants for the analysis.

sumFun<-function(dat,vars,type,digs=2,na.rm=TRUE){
  # dat = data frame to be summarised
  # vars = character vector with variables names
  # type = character vector of same length as vars specifying the type of summary to compute (needs to be one of "mean_sd","median_iqr","n_perc")
  # digs = integer specifying the number of decimal digits to report (defaults to 2)
  # na.rm = whether to remove NAs or not (defaults to TRUE)
  
  res<-data.frame(
    var=character(0),
    stat1=character(0),
    stat2=character(0)
  )
  
  for(i in 1:length(vars)){
    if(type[i]=="mean_sd"){
      stat1<-format(nsmall=digs,round(digits=digs,mean(dat[,vars[i]],na.rm=T)))
      stat2<-paste(sep="","(",format(nsmall=digs,round(digits=digs,sd(dat[,vars[i]],na.rm=T))),")")
      #varName<-paste(sep=" ",vars[i],"; mean (sd)")
      varName<-paste(sep=" ","mean (sd)")
      
      res<-rbind(res,c(varName,stat1,stat2))
      rm(varName,stat1,stat2)
    }else if(type[i]=="median_iqr"){
      stat1<-format(nsmall=digs,round(digits=digs,median(dat[,vars[i]],na.rm=T)))
      stat2<-paste(sep="","(",paste(collapse=",",format(nsmall=digs,round(digits=digs,quantile(dat[,vars[i]],probs=c(0.25,0.75),na.rm=T)))),")")
      
      #varName<-paste(sep=" ",vars[i],"; median (IQR)")
      varName<-paste(sep=" ","median (IQR)")
      res<-rbind(res,c(varName,stat1,stat2))
      rm(varName,stat1,stat2)
    }else if(type[i]=="n_perc"){
      for(val in sort(unique(dat[,vars[i]]))){
        stat1<-sum(dat[,vars[i]]==val,na.rm=T)
        stat2<-paste(sep="","(",format(nsmall=digs,round(digits=digs,100*sum(dat[,vars[i]]==val,na.rm=T)/sum(!is.na(dat[,vars[i]])))),"%)")
        
        #varName<-paste(sep="",vars[i],"; ",val,"; n (%)")
        varName<-paste(sep="",val,"; n (%)")
        res<-rbind(res,c(varName,stat1,stat2))
        rm(varName,stat1,stat2)
      }
    }
  }
  
  colnames(res)<-c("variable","stat1","stat2")
  return(res)
}
  
datSum<-sumFun(dat,
               vars=c("hosp_com","site","sex","age","bmi","te","apri","gpr","fib4","alt"),
               type=c("n_perc","n_perc","n_perc","median_iqr","mean_sd","median_iqr","median_iqr","median_iqr","median_iqr","median_iqr","median_iqr"))

datSum %>%
  kable(caption="Patient characteristics",row.names=F) %>%
  kable_styling(full_width=F) %>%
  pack_rows(start_row=1,end_row=2,group_label=paste(sep="","Type of study; n = ",sum(!is.na(dat$hosp_com))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$hosp_com))/nrow(dat))),"%)")) %>%
  pack_rows(start_row=3,end_row=13,group_label=paste(sep="","Site; n = ",sum(!is.na(dat$site))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$site))/nrow(dat))),"%)")) %>%
  pack_rows(start_row=14,end_row=15,group_label=paste(sep="","Sex; n = ",sum(!is.na(dat$sex))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$age))/nrow(dat))),"%)")) %>%
  pack_rows(start_row=16,end_row=16,group_label=paste(sep="","Age; n = ",sum(!is.na(dat$age))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$age))/nrow(dat))),"%)")) %>%
  pack_rows(start_row=17,end_row=17,group_label=paste(sep="","BMI; n = ",sum(!is.na(dat$bmi))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$bmi))/nrow(dat))),"%)")) %>%
  pack_rows(start_row=18,end_row=18,group_label=paste(sep="","Transient elastography; n = ",sum(!is.na(dat$te))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$te))/nrow(dat))),"%)")) %>%
  pack_rows(start_row=19,end_row=19,group_label=paste(sep="","APRI (AST to platelet ratio); n = ",sum(!is.na(dat$apri))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$apri))/nrow(dat))),"%)")) %>%
  pack_rows(start_row=20,end_row=20,group_label=paste(sep="","GPR (GGT to platelet ratio); n = ",sum(!is.na(dat$gpr))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$gpr))/nrow(dat))),"%)")) %>%
  pack_rows(start_row=21,end_row=21,group_label=paste(sep="","Fibrosis 4 score; n = ",sum(!is.na(dat$bmi))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$bmi))/nrow(dat))),"%)")) %>%
  pack_rows(start_row=22,end_row=22,group_label=paste(sep="","ALT; n = ",sum(!is.na(dat$bmi))," (",format(nsmall=1,round(digits=1,100*sum(!is.na(dat$bmi))/nrow(dat))),"%)"))
Table 1: Patient characteristics
variable stat1 stat2
Type of study; n = 2945 (100.0%)
Community; n (%) 735 (24.96%)
Hospital; n (%) 2210 (75.04%)
Site; n = 2945 (100.0%)
Burkina Faso; n (%) 35 (1.19%)
Cape Town; n (%) 155 (5.26%)
Ethiopia; n (%) 1047 (35.55%)
Gambia; n (%) 782 (26.55%)
Malawi; n (%) 100 (3.40%)
Nigeria; n (%) 190 (6.45%)
Senegal_HPB; n (%) 84 (2.85%)
Senegal_Maud; n (%) 275 (9.34%)
Senegal_PRF; n (%) 151 (5.13%)
Stellenbosch; n (%) 36 (1.22%)
Zambia; n (%) 90 (3.06%)
Sex; n = 2945 (100.0%)
Female; n (%) 1138 (38.64%)
Male; n (%) 1807 (61.36%)
Age; n = 2945 (100.0%)
median (IQR) 34.00 (28.00,41.00)
BMI; n = 2583 (87.7%)
mean (sd) 23.12 (4.71)
Transient elastography; n = 2945 (100.0%)
median (IQR) 5.60 (4.50,7.10)
APRI (AST to platelet ratio); n = 2945 (100.0%)
median (IQR) 0.32 (0.22,0.48)
GPR (GGT to platelet ratio); n = 2261 (76.8%)
median (IQR) 0.12 (0.07,0.18)
Fibrosis 4 score; n = 2583 (87.7%)
median (IQR) 0.85 (0.58,1.30)
ALT; n = 2583 (87.7%)
median (IQR) 25.00 (19.00,35.00)

Definiton of the gold standard diagnosis

Reference standards (which are assumed to be gold standards) are based on variable te, transient elastrography result. Specifically:

In the dataset, these 3 gold standards are encoded by variables te_115, te_95 and te_79.

dat<-dat %>%
  mutate(
    te_115=ifelse(te>11.5,1,0),
    te_95=ifelse(te>9.5,1,0),
    te_79=ifelse(te>7.9,1,0),
  )

Diagnostics to be evaluated

dat<-dat %>%
  mutate(
    apri_20=ifelse(apri>=2,1,0),
    apri_15=ifelse(apri>=1.5,1,0),
    gpr_56=ifelse(gpr>=0.56,1,0),
    gpr_32=ifelse(gpr>=0.32,1,0),
  )

Very naive analysis (to get feel for data)

Let’s just assume we could analyse all individual-level data as is, simply pulling out the sensitivity, specificity rather than fitting a model that considers random effects and individual- and study-level covariates.

gr<-rbind(
  c(11.5,"2.0","0.56"),
  c(9.5,"2.0","0.56"),
  c(7.9,"1.5","0.32")
)

rnames<-character(0)
for(s in c("All",sort(unique(dat$Sitecountry)))){
  for(i in 1:3){
    rnames<-c(rnames,paste(sep="_",paste(sep="","Site",s),paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))))
  }
}
rnames<-gsub(rnames,pattern="SiteAll",replacement="All")

sensSpecMat<-matrix(nrow=3*(length(unique(dat$Sitecountry))+1),ncol=8)
rownames(sensSpecMat)<-rnames
colnames(sensSpecMat)<-c("site","te","apri","gpr","apri_sens","apri_spec","gpr_sens","gpr_spec")
sensSpecMat<-as.data.frame(sensSpecMat)

for(s in c("All",sort(unique(dat$Sitecountry)))){
  if(s=="All"){
    datTmp<-dat
  }else{
    datTmp<-dat %>% filter(Sitecountry==s)
    s<-paste(sep="","Site",s)
  }
    
  for(i in 1:nrow(gr)){
    gthr<-gr[i,1]
    thr1<-gr[i,2]
    thr2<-gr[i,3]
    
    sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"site"]<-gsub(s,pattern="Site",replacement="")
    sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"te"]<-gthr
    sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"apri"]<-thr1
    sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"gpr"]<-thr2
    
    sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"apri_sens"]<-mean(datTmp[datTmp[,paste(sep="_","te",gsub(gthr,pattern="\\.",replacement=""))]==1,paste(sep="_","apri",gsub(thr1,pattern="\\.",replacement=""))])
    sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"gpr_sens"]<-mean(datTmp[datTmp[,paste(sep="_","te",gsub(gthr,pattern="\\.",replacement=""))]==1,paste(sep="_","gpr",gsub(thr2,pattern="0\\.",replacement=""))],na.rm=T)
    
    sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"apri_spec"]<-mean(1-datTmp[datTmp[,paste(sep="_","te",gsub(gthr,pattern="\\.",replacement=""))]==0,paste(sep="_","apri",gsub(thr1,pattern="\\.",replacement=""))])
    sensSpecMat[paste(sep="_",s,paste(collapse="_",paste(sep="",c("te","apri","gpr"),gr[i,]))),"gpr_spec"]<-mean(1-datTmp[datTmp[,paste(sep="_","te",gsub(gthr,pattern="\\.",replacement=""))]==0,paste(sep="_","gpr",gsub(thr2,pattern="0\\.",replacement=""))],na.rm=T)
    
  }
}

sensSpecMatKable<-sensSpecMat
for(i in 1:nrow(sensSpecMat)){
  for(j in 5:ncol(sensSpecMat)){
    if(!is.na(sensSpecMat[i,j])){
      sensSpecMatKable[i,j]<-paste(sep="",format(nsmall=1,round(digits=1,100*sensSpecMat[i,j])),"%")
    }else{
      sensSpecMatKable[i,j]<-"-"
    }
  }
}
rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te11.5_apri2.0_gpr0.56")]<-gsub(rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te11.5_apri2.0_gpr0.56")],pattern="_te11.5_apri2.0_gpr0.56",replacement="; te > 11.5 (apri > 2.0, gpr > 0.56)")
rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te9.5_apri2.0_gpr0.56")]<-gsub(rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te9.5_apri2.0_gpr0.56")],pattern="_te9.5_apri2.0_gpr0.56",replacement="; te > 9.5 (apri > 2.0, gpr > 0.56)")
rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te7.9_apri1.5_gpr0.32")]<-gsub(rownames(sensSpecMatKable)[grepl(rownames(sensSpecMatKable),pattern="_te7.9_apri1.5_gpr0.32")],pattern="_te7.9_apri1.5_gpr0.32",replacement="; te > 7.9 (apri > 1.5, gpr > 0.32)")

sensSpecMatKable <- cbind("",sensSpecMatKable)

sensSpecMatKable %>%
  dplyr::select(!site) %>%
  kable(caption="Naive sensitivity and specificity estimates (expressed as percentages).",col.names=c(" ","TE","APRI","GPR","Sensitivity","Specificity","Sensitivity","Specificity"),row.names = FALSE) %>%
  kable_styling(full_width = FALSE) %>%
  add_header_above(c(" "=1,"Diagnostic thresholds"=3,"APRI"=2,"GPR"=2)) %>%
  pack_rows(start_row=1,end_row=3,group_label="All samples") %>%
  pack_rows(start_row=4,end_row=6,group_label="Site 1") %>%
  pack_rows(start_row=7,end_row=9,group_label="Site 2") %>%
  pack_rows(start_row=10,end_row=12,group_label="Site 3") %>%
  pack_rows(start_row=13,end_row=15,group_label="Site 4") %>%
  pack_rows(start_row=16,end_row=18,group_label="Site 5") %>%
  pack_rows(start_row=19,end_row=21,group_label="Site 6") %>%
  pack_rows(start_row=22,end_row=24,group_label="Site 7") %>%
  pack_rows(start_row=25,end_row=27,group_label="Site 8") %>%
  pack_rows(start_row=28,end_row=30,group_label="Site 9") %>%
  pack_rows(start_row=31,end_row=33,group_label="Site 10") %>%
  pack_rows(start_row=34,end_row=36,group_label="Site 11")
Table 2: Naive sensitivity and specificity estimates (expressed as percentages).
Diagnostic thresholds
APRI
GPR
TE APRI GPR Sensitivity Specificity Sensitivity Specificity
All samples
11.5 2.0 0.56 20.2% 99.5% 31.6% 98.0%
9.5 2.0 0.56 15.5% 99.5% 26.5% 98.3%
7.9 1.5 0.32 14.4% 99.2% 37.1% 94.6%
Site 1
11.5 2.0 0.56 10.5% 99.8% 18.6% 99.7%
9.5 2.0 0.56 8.5% 99.8% 15.1% 99.7%
7.9 1.5 0.32 10.1% 99.3% 29.3% 97.5%
Site 2
11.5 2.0 0.56 0.0% 99.3% 22.2% 95.0%
9.5 2.0 0.56 0.0% 99.3% 20.0% 96.2%
7.9 1.5 0.32 5.3% 99.1% 44.7% 92.9%
Site 3
11.5 2.0 0.56
100.0%
100.0%
9.5 2.0 0.56 0.0% 100.0% 0.0% 100.0%
7.9 1.5 0.32 0.0% 100.0% 20.0% 98.2%
Site 4
11.5 2.0 0.56 70.0% 100.0% 83.3% 98.6%
9.5 2.0 0.56 65.6% 100.0% 78.1% 98.5%
7.9 1.5 0.32 63.2% 100.0% 84.2% 95.2%
Site 5
11.5 2.0 0.56 26.1% 100.0%
9.5 2.0 0.56 18.2% 100.0%
7.9 1.5 0.32 18.9% 99.3%
Site 6
11.5 2.0 0.56 40.0% 100.0%
9.5 2.0 0.56 28.6% 100.0%
7.9 1.5 0.32 23.1% 100.0%
Site 7
11.5 2.0 0.56 33.3% 97.5% 66.7% 86.1%
9.5 2.0 0.56 26.7% 98.7% 66.7% 90.9%
7.9 1.5 0.32 25.0% 98.5% 44.4% 76.7%
Site 8
11.5 2.0 0.56 21.4% 99.1% 42.9% 97.4%
9.5 2.0 0.56 10.0% 99.1% 27.6% 97.6%
7.9 1.5 0.32 4.1% 98.6% 31.4% 91.6%
Site 9
11.5 2.0 0.56 13.3% 100.0% 33.3% 98.5%
9.5 2.0 0.56 11.8% 100.0% 29.4% 98.5%
7.9 1.5 0.32 10.0% 100.0% 50.0% 95.4%
Site 10
11.5 2.0 0.56 100.0% 100.0% 100.0% 100.0%
9.5 2.0 0.56 50.0% 100.0% 50.0% 100.0%
7.9 1.5 0.32 50.0% 100.0% 50.0% 100.0%
Site 11
11.5 2.0 0.56 8.7% 99.2% 9.1% 97.5%
9.5 2.0 0.56 9.4% 99.6% 13.3% 98.3%
7.9 1.5 0.32 10.2% 99.5% 25.0% 97.1%

Question:

Sensitivity and specificity are extremely variable from study to study. I find it difficult to understand this. You would expect some heterogeneity due to how certain diagnostics are delivered, some underlying characteristics that impct test results which vary from population to population, but not quite to this extent. Positive and negative predictive values can vary a lot from one population to another, but that is because these, unlike sensitvity and specificity, depend on the prevalence of the condition under investigation.

Would be good to discuss this.

Model

Earlier versions of this document (20210627 and earlier) compared several different models. From this, the decision was made for a model with individual- and study-level covariates and study-specific random effects

This is model (4) from Riley et al (2008).

This model estimates summary sensitivity and specificity with study random effects and individual- and study-level covariates.

We include the following co-variates:

Pre-defined thresholds for APRI and GPR

As noted earlier in this document, there exist recommended thresholds for APRI (>2.0 for F4 and >1.5 for F2) and GPR (>0.56 for F4 and >0.32 for F2). Together with the 2 thresholds for the reference standard for F4 (te>11.5 and te>9.5) and 1 threshold for F2 (te>7.9), this gives 6 models for pre-defined models. The overall sensitivities and specificities for these models ar reported in Tables 3 and 4 below. Model-specific tables giving ORs for the different co-variates are given in the appendix.

datJags<-list()

datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_115==1,1,2)
datJags$y<-dat$apri_20
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc

randSeed<-20210615
set.seed(randSeed)

jagsModel4_te115_apri20 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te115_apri20<-coda.samples(model=jagsModel4_te115_apri20,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)

tmp<-try(summary(parsModel4_te115_apri20)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te115_apri20))
mod_te115_apri20<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()

datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_95==1,1,2)
datJags$y<-dat$apri_20
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc

randSeed<-20210615
set.seed(randSeed)

jagsModel4_te95_apri20 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te95_apri20<-coda.samples(model=jagsModel4_te95_apri20,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)

tmp<-try(summary(parsModel4_te95_apri20)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te95_apri20))
mod_te95_apri20<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()

datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_79==1,1,2)
datJags$y<-dat$apri_15
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc

randSeed<-20210615
set.seed(randSeed)

jagsModel4_te79_apri15 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te79_apri15<-coda.samples(model=jagsModel4_te79_apri15,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)

tmp<-try(summary(parsModel4_te79_apri15)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te79_apri15))
mod_te79_apri15<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()

datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_115==1,1,2)
datJags$y<-dat$gpr_56
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc

randSeed<-20210615
set.seed(randSeed)

jagsModel4_te115_gpr56 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te115_gpr56<-coda.samples(model=jagsModel4_te115_gpr56,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)

tmp<-try(summary(parsModel4_te115_gpr56)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te115_gpr56))
mod_te115_gpr56<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()

datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_95==1,1,2)
datJags$y<-dat$gpr_56
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc

randSeed<-20210615
set.seed(randSeed)

jagsModel4_te95_gpr56 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 10000) # needs to be large!
parsModel4_te95_gpr56<-coda.samples(model=jagsModel4_te95_gpr56,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)

tmp<-try(summary(parsModel4_te95_gpr56)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te95_gpr56))
mod_te95_gpr56<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)
datJags<-list()

datJags$N<-nrow(dat)
datJags$m<-max(as.integer(factor(dat$Sitecountry)))
datJags$study<-as.integer(factor(dat$Sitecountry))
datJags$state<-ifelse(dat$te_79==1,1,2)
datJags$y<-dat$gpr_32
datJags$studyTypeHospital<-dat$studyTypeHospital
datJags$testReasonL<-dat$testReasonL
datJags$alc<-dat$alc

randSeed<-20210615
set.seed(randSeed)

jagsModel4_te79_gpr32 <- jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
parsModel4_te79_gpr32<-coda.samples(model=jagsModel4_te79_gpr32,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)

tmp<-try(summary(parsModel4_te79_gpr32)$statistics,silent=TRUE)
tmpCI<-t(hdi(parsModel4_te79_gpr32))
mod_te79_gpr32<-cbind(tmp,tmpCI)
rm(tmp,tmpCI)

save(list=c("jagsModel4_te115_apri20","jagsModel4_te95_apri20","jagsModel4_te79_apri15","jagsModel4_te115_gpr56","jagsModel4_te95_gpr56","jagsModel4_te79_gpr32","parsModel4_te115_apri20","parsModel4_te95_apri20","parsModel4_te79_apri15","parsModel4_te115_gpr56","parsModel4_te95_gpr56","parsModel4_te79_gpr32"),file=paste(sep="","fixedThresholdModelsResultsList_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
fixedThrSumTabApri %>%
  knitr::kable(caption = "Summaries of the sensitivites and specificities for the 3 models for pre-defined thresholds on APRI.", col.names=c("te>11.5","te>9.5","te>7.9"),row.names = FALSE,align = "c") %>%
  add_header_above(c("APRI>=2.0"=1,"APRI>=2.0"=1,"APRI>=1.5"=1)) %>%
  pack_rows(group_label="sensitivity\n(95% CI)",1,2) %>%
  pack_rows(group_label="specificity\n(95% CI)",3,4) %>%
  kable_styling(full_width = FALSE)
Table 3: Summaries of the sensitivites and specificities for the 3 models for pre-defined thresholds on APRI.
APRI>=2.0
APRI>=2.0
APRI>=1.5
te>11.5 te>9.5 te>7.9
sensitivity
(95% CI)
28.76% 14.76% 7.32%
(5.70%,54.83%) (2.70%,28.52%) (1.87%,13.93%)
specificity
(95% CI)
99.18% 99.16% 99.16%
(98.20%,99.95%) (98.18%,99.90%) (98.17%,99.93%)
Table 4: Summaries of the sensitivites and specificities for the 3 models for pre-defined thresholds on GPR.
GPR>=0.56
GPR>=0.56
GPR>=0.32
te>11.5 te>9.5 te>7.9
sensitivity
(95% CI)
22.68% 20.43% 18.88%
(3.33%,47.31%) (5.34%,37.95%) (8.12%,30.73%)
specificity
(95% CI)
98.45% 98.88% 96.76%
(97.11%,99.56%) (97.88%,99.69%) (94.83%,98.50%)

Optimum thresholds for APRI

bestModels<-list()
bestThrs<-list()
aucs<-list()
doForest<-function(dat,outFile,marker,teThr,bestThr){
  datSumForest<-dat %>%
    group_by(site) %>%
    filter(!is.na(get(marker))) %>%
    summarise(
      n=n(),
      nWithCompleteData=sum(!is.na(te) & !is.na(get(marker)) & !is.na(testReasonL) & !is.na(studyTypeHospital) & !is.na(alc)),
      nTePos=sum(te>teThr),
      propTePos=format(nsmall=4,round(digits=4,mean(te>teThr))),
      nMarkerPos=sum(get(marker)>bestThr),
      TP=sum(te>teThr & get(marker)>=bestThr),
      FN=sum(te>teThr & get(marker)<bestThr),
      FP=sum(te<=teThr & get(marker)>=bestThr),
      TN=sum(te<=teThr & get(marker)<bestThr)
    ) %>%
    mutate(
      sens=TP/(TP+FN),
      sensCI=NA,
      sensLower=NA,
      sensUpper=NA,
      spec=TN/(FP+TN),
      specCI=NA,
      specLower=NA,
      specUpper=NA,
      issum=FALSE
    )
  
  for(j in 1:nrow(datSumForest)){
    if(datSumForest$TP[j]+datSumForest$FN[j]>0){
      datSumForest$sensCI[j]=paste(sep="","(",paste(collapse=",",format(nsmall=2,round(digits=2,binom.test(x=datSumForest$TP[j],n=datSumForest$TP[j]+datSumForest$FN[j])$conf.int))),")")
      datSumForest$sensLower[j]=binom.test(x=datSumForest$TP[j],n=datSumForest$TP[j]+datSumForest$FN[j])$conf.int[1]
      datSumForest$sensUpper[j]=binom.test(x=datSumForest$TP[j],n=datSumForest$TP[j]+datSumForest$FN[j])$conf.int[2]
    }
    if(datSumForest$FP[j]+datSumForest$TN[j]>0){
      datSumForest$specCI[j]=paste(sep="","(",paste(collapse=",",format(nsmall=2,round(digits=2,binom.test(x=datSumForest$TN[j],n=datSumForest$FP[j]+datSumForest$TN[j])$conf.int))),")")
      datSumForest$specLower[j]=binom.test(x=datSumForest$TN[j],n=datSumForest$FP[j]+datSumForest$TN[j])$conf.int[1]
      datSumForest$specUpper[j]=binom.test(x=datSumForest$TN[j],n=datSumForest$FP[j]+datSumForest$TN[j])$conf.int[2]
    }
  }
  
  sumForest <- data.frame(
    site="Summary",
    n=sum(datSumForest$n), 
    nWithCompleteData=sum(datSumForest$nWithCompleteData),
    nTePos=sum(datSumForest$nTePos),
    propTePos=format(nsmall=4,round(digits=4,sum(datSumForest$nTePos)/sum(datSumForest$n))),
    nMarkerPos=sum(datSumForest$nMarkerPos),
    TP=sum(datSumForest$TP),
    FN=sum(datSumForest$FN),
    FP=sum(datSumForest$FP),
    TN=sum(datSumForest$TN),
    sens = round(digits=4,bestModels[[paste(sep="",marker,"_te",10*teThr)]]["overall sensitivity","Mean"]),
    sensCI = bestModels[[paste(sep="",marker,"_te",10*teThr)]]["overall sensitivity","CrI"],
    sensLower=as.numeric(gsub(pattern="\\(",replacement="",unlist(strsplit(split=",",bestModels[[paste(sep="",marker,"_te",10*teThr)]]["overall sensitivity","CrI"]))[1])),
    sensUpper=as.numeric(gsub(pattern="\\)",replacement="",unlist(strsplit(split=",",bestModels[[paste(sep="",marker,"_te",10*teThr)]]["overall sensitivity","CrI"]))[2])),
    spec = round(digits=4,bestModels[[paste(sep="",marker,"_te",10*teThr)]]["overall specificity","Mean"]),
    specCI = bestModels[[paste(sep="",marker,"_te",10*teThr)]]["overall specificity","CrI"],
    specLower=as.numeric(gsub(pattern="\\(",replacement="",unlist(strsplit(split=",",bestModels[[paste(sep="",marker,"_te",10*teThr)]]["overall specificity","CrI"]))[1])),
    specUpper=as.numeric(gsub(pattern="\\)",replacement="",unlist(strsplit(split=",",bestModels[[paste(sep="",marker,"_te",10*teThr)]]["overall specificity","CrI"]))[2])),
    issum = TRUE)
  
  sumForestHeader <- data.frame(
    site = c("","Site"),
    n = c("n","(all)"),
    nWithCompleteData=c("n","(complete data)"),
    nTePos=c("n",paste(sep="","(te > ",teThr,")")),
    propTePos=c("prop",paste(sep="","(te > ",teThr,")")),
    nMarkerPos=c("n",paste(sep="","(",toupper(marker),">=",bestThr,")")),
    TP=c("","TP"),
    FN=c("","FN"),
    FP=c("","FP"),
    TN=c("","TN"),
    sens=c("sensitivity",""),
    sensCI = c("sensitivity","(95% CI)"),
    spec=c("specificity",""),
    specCI = c("specificity","(95% CI)"),
    issum = TRUE)
  
  sumForestEmpty <- data.frame(mean = NA_real_)
  
  datSumForest<-datSumForest %>%
    mutate(
      n=as.character(n),
      nWithCompleteData=as.character(nWithCompleteData),
      nTePos=as.character(nTePos),
      propTePos=as.character(propTePos),
      nMarkerPos=as.character(nMarkerPos),
      TP=as.character(TP),
      FN=as.character(FN),
      FP=as.character(FP),
      TN=as.character(TN),
      sens=format(nsmall=2,round(digits=2,(sens))),
      spec=format(nsmall=2,round(digits=2,(spec)))
    )
  
  sumForest<-sumForest %>%
    mutate(
      n=as.character(n),
      nWithCompleteData=as.character(nWithCompleteData),
      nTePos=as.character(nTePos),
      propTePos=as.character(propTePos),
      nMarkerPos=as.character(nMarkerPos),
      TP=as.character(TP),
      FN=as.character(FN),
      FP=as.character(FP),
      TN=as.character(TN),
      sens=format(nsmall=4,round(digits=4,(sens))),
      spec=format(nsmall=4,round(digits=4,(spec)))
    )
  
  forestOut <- bind_rows(sumForestHeader,
                         datSumForest,
                         #sumForestEmpty,
                         sumForest)
  
  forestOut$sens[is.nan(forestOut$sens) | grepl(forestOut$sens,pattern="NaN")]<-NA
  
  gSens<-grid.grabExpr(print(forestOut %>%
    mutate(
      mean=as.numeric(sens),
      lower=sensLower,
      upper=sensUpper
    ) %>%
    dplyr::select(site,n,nTePos,nMarkerPos,TP,FN,FP,TN,sens,sensCI,mean,lower,upper) %>%
    as_tibble() %>%
    forestplot(
      labeltext=c(site,n,nTePos,nMarkerPos,TP,FN,FP,TN,sens,sensCI),
      is.summary = forestOut$issum,
      col = fpColors(box = "royalblue",
                     line = "darkblue",
                     summary = "royalblue"),
      graphwidth=unit(6,"cm"),
      clip=c(0,1))
  ))  
    
  gSpec<-grid.grabExpr(print(forestOut %>%
    mutate(
      mean=as.numeric(spec),
      lower=specLower,
      upper=specUpper
    ) %>%
    dplyr::select(spec,specCI,mean,lower,upper) %>%
    as_tibble() %>%
    forestplot(
      labeltext=c(spec,specCI),
      is.summary = forestOut$issum,
      col = fpColors(box = "royalblue",
                     line = "darkblue",
                     summary = "royalblue"),
      graphwidth=unit(6,"cm"),
      clip=c(0,1))
  ))
  
  png(width=19,height=8,units="in",res=300,file=paste(sep="","forestPlot",str_to_title(marker),"Te",10*teThr,gsub(Sys.Date(),pattern="-",replacement=""),".png"))
  grid.arrange(gSens,gSpec,nrow=1,widths=c(13.5,5.5))
  dev.off()
}

F4 - te > 11.5kPa

apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43

apri115ResultsList<-list()
apri115ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(apriVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  apr<-apriVect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_115==1,1,2)
  datJags$y<-ifelse(dat$apri>=apr,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(apriVect)){
  apr<-apriVect[i]
  
  apri115ModList[[paste(sep="","apri_mod_",apr)]]<-parList[[i]]$model
  apri115ModList[[paste(sep="","apri_samples_",apr)]]<-parList[[i]]$pars
  apri115ResultsList[[paste(sep="","apri",apr)]]<-parList[[i]]$results
}

rm(parList)
save(apri115ResultsList,apri115ModList,file=paste(sep="","apriResultsList_te115_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43

load("apriResultsList_te115_20210805.RData")
sensVect<-numeric(length(apriVect))
specVect<-numeric(length(apriVect))
for(i in 1:length(apriVect)){
  sensVect[i]<-apri115ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-apri115ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  apri=apriVect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-apriVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nAPRI threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for APRI with diagnosis threshold te > 11.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on APRI (from 0 to 13.26 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

apri115ResultsList[[paste(sep="","apri",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) APRI model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 5: Best (according to Youden’s J) APRI model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 0.5804 (0.0958,1.2771)
hazardous alcohol drinking (OR for specificity) 1.0163 (0.5602,1.5236)
suspected liver disease screening (OR for sensitivity) 3.7027 (0.6861,8.2046)
suspected liver disease screening (OR for specificity) 0.1457 (0.0788,0.2173)
hospital-based study (OR for sensitivity) 2.4776 (0.2120,6.5763)
hospital-based study (OR for specificity) 1.4508 (1.0241,1.9267)
random effects variance (logit sensitivity) 1.2545 (0.1256,3.4243)
random effects variance (logit specificity) 0.7847 (0.2034,1.6496)
overall sensitivity 0.6222 (0.3468,0.8748)
overall specificity 0.8564 (0.8162,0.8962) 
bestModels[["apri_te115"]]<-apri115ResultsList[[paste(sep="","apri",bestThr)]]
bestThrs[["apri_te115"]]<-bestThr
aucs[["apri_te115"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="apri",teThr=11.5)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotApriTe115",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for APRI with te >= 11.5kPa.

Figure 1: Forest plots for APRI with te >= 11.5kPa.

F4 - te > 9.5kPa

apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43

apri95ResultsList<-list()
apri95ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(apriVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  apr<-apriVect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_95==1,1,2)
  datJags$y<-ifelse(dat$apri>=apr,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(apriVect)){
  apr<-apriVect[i]
  
  apri95ModList[[paste(sep="","apri_mod_",apr)]]<-parList[[i]]$model
  apri95ModList[[paste(sep="","apri_samples_",apr)]]<-parList[[i]]$pars
  apri95ResultsList[[paste(sep="","apri",apr)]]<-parList[[i]]$results
}

rm(parList)
save(apri95ResultsList,apri95ModList,file=paste(sep="","apriResultsList_te95_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43

load("apriResultsList_te95_20210805.RData")
sensVect<-numeric(length(apriVect))
specVect<-numeric(length(apriVect))
for(i in 1:length(apriVect)){
  sensVect[i]<-apri95ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-apri95ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  apri=apriVect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-apriVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nAPRI threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for APRI with diagnosis threshold te > 9.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on APRI (from 0 to 13.26 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

apri95ResultsList[[paste(sep="","apri",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) APRI model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 6: Best (according to Youden’s J) APRI model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 0.3891 (0.0839, 0.8028)
hazardous alcohol drinking (OR for specificity) 1.0314 (0.5904, 1.5065)
suspected liver disease screening (OR for sensitivity) 13.3700 (3.6185,26.6388)
suspected liver disease screening (OR for specificity) 0.1644 (0.1054, 0.2328)
hospital-based study (OR for sensitivity) 0.9215 (0.2037, 1.8660)
hospital-based study (OR for specificity) 1.4098 (1.0479, 1.7939)
random effects variance (logit sensitivity) 0.7792 (0.1152, 1.9096)
random effects variance (logit specificity) 0.7465 (0.1985, 1.5735)
overall sensitivity 0.7049 (0.5192, 0.8895)
overall specificity 0.7750 (0.7263, 0.8235) 
bestModels[["apri_te95"]]<-apri95ResultsList[[paste(sep="","apri",bestThr)]]
bestThrs[["apri_te95"]]<-bestThr
aucs[["apri_te95"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="apri",teThr=9.5)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotApriTe95",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for APRI with te >= 9.5kPa.

Figure 2: Forest plots for APRI with te >= 9.5kPa.

F2 - te > 7.9kPa

apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43

apri79ResultsList<-list()
apri79ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(apriVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  apr<-apriVect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_79==1,1,2)
  datJags$y<-ifelse(dat$apri>=apr,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(apriVect)){
  apr<-apriVect[i]
  
  apri79ModList[[paste(sep="","apri_mod_",apr)]]<-parList[[i]]$model
  apri79ModList[[paste(sep="","apri_samples_",apr)]]<-parList[[i]]$pars
  apri79ResultsList[[paste(sep="","apri",apr)]]<-parList[[i]]$results
}

rm(parList)
save(apri79ResultsList,apri79ModList,file=paste(sep="","apriResultsList_te79_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
apriVect<-unique(c(0,1.5,2,round(digits=2,quantile(dat$apri,probs=seq(0,1,length=40)))))
apriVect<-sort(apriVect) # length 43

load("apriResultsList_te79_20210805.RData")
sensVect<-numeric(length(apriVect))
specVect<-numeric(length(apriVect))
for(i in 1:length(apriVect)){
  sensVect[i]<-apri79ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-apri79ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  apri=apriVect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-apriVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nAPRI threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for APRI with diagnosis threshold te > 7.9kPa (F2).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on APRI (from 0 to 13.26 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

apri79ResultsList[[paste(sep="","apri",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) APRI model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 7: Best (according to Youden’s J) APRI model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 0.3985 (0.1254, 0.7339)
hazardous alcohol drinking (OR for specificity) 0.9596 (0.5289, 1.4377)
suspected liver disease screening (OR for sensitivity) 10.7007 (4.5094,18.3832)
suspected liver disease screening (OR for specificity) 0.1582 (0.0914, 0.2302)
hospital-based study (OR for sensitivity) 0.6865 (0.2929, 1.1543)
hospital-based study (OR for specificity) 1.3991 (1.0028, 1.8296)
random effects variance (logit sensitivity) 0.5718 (0.1087, 1.3040)
random effects variance (logit specificity) 0.9175 (0.2338, 1.9148)
overall sensitivity 0.6674 (0.5289, 0.8061)
overall specificity 0.8091 (0.7597, 0.8563) 
bestModels[["apri_te79"]]<-apri79ResultsList[[paste(sep="","apri",bestThr)]]
bestThrs[["apri_te79"]]<-bestThr
aucs[["apri_te79"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="apri",teThr=7.9)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotApriTe79",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for APRI with te >= 7.9kPa.

Figure 3: Forest plots for APRI with te >= 7.9kPa.

Optimum thresholds for GPR

Note that the GPR data is the only one of APRI, GPR, ALT, FIB4 to have missing values. These are ignored when selecting the threshold points. The Bayesian models will integrate them out in a principled fashion.

F4 - te > 11.5kPa

gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40

gpr115ResultsList<-list()
gpr115ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(gprVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  gpr<-gprVect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_115==1,1,2)
  datJags$y<-ifelse(dat$gpr>=gpr,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(gprVect)){
  gpr<-gprVect[i]
  
  gpr115ModList[[paste(sep="","gpr_mod_",gpr)]]<-parList[[i]]$model
  gpr115ModList[[paste(sep="","gpr_samples_",gpr)]]<-parList[[i]]$pars
  gpr115ResultsList[[paste(sep="","gpr",gpr)]]<-parList[[i]]$results
}

rm(parList)
save(gpr115ResultsList,gpr115ModList,file=paste(sep="","gprResultsList_te115_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40

load("gprResultsList_te115_20210805.RData")
sensVect<-numeric(length(gprVect))
specVect<-numeric(length(gprVect))
for(i in 1:length(gprVect)){
  sensVect[i]<-gpr115ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-gpr115ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  gpr=gprVect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-gprVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nGPR threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for GPR with diagnosis threshold te > 11.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on GPR (from 0 to 13.65 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

gpr115ResultsList[[paste(sep="","gpr",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) GPR model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 8: Best (according to Youden’s J) GPR model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 1.1523 (0.0577, 3.2957)
hazardous alcohol drinking (OR for specificity) 0.5383 (0.3125, 0.8002)
suspected liver disease screening (OR for sensitivity) 11.3094 (0.3977,34.4099)
suspected liver disease screening (OR for specificity) 0.3301 (0.2145, 0.4602)
hospital-based study (OR for sensitivity) 1.8280 (0.0087, 5.5713)
hospital-based study (OR for specificity) 0.4535 (0.3272, 0.5837)
random effects variance (logit sensitivity) 1.1528 (0.0905, 3.3541)
random effects variance (logit specificity) 0.4853 (0.1088, 1.0596)
overall sensitivity 0.8622 (0.6569, 0.9994)
overall specificity 0.6979 (0.6299, 0.7637) 
bestModels[["gpr_te115"]]<-gpr115ResultsList[[paste(sep="","gpr",bestThr)]]
bestThrs[["gpr_te115"]]<-bestThr
aucs[["gpr_te115"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="gpr",teThr=11.5)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotGprTe115",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for GPR with te >= 11.5kPa.

Figure 4: Forest plots for GPR with te >= 11.5kPa.

F4 - te > 9.5kPa

gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40

gpr95ResultsList<-list()
gpr95ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(gprVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  gpr<-gprVect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_95==1,1,2)
  datJags$y<-ifelse(dat$gpr>=gpr,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(gprVect)){
  gpr<-gprVect[i]
  
  gpr95ModList[[paste(sep="","gpr_mod_",gpr)]]<-parList[[i]]$model
  gpr95ModList[[paste(sep="","gpr_samples_",gpr)]]<-parList[[i]]$pars
  gpr95ResultsList[[paste(sep="","gpr",gpr)]]<-parList[[i]]$results
}

rm(parList)
save(gpr95ResultsList,gpr95ModList,file=paste(sep="","gprResultsList_te95_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40

load("gprResultsList_te95_20210805.RData")
sensVect<-numeric(length(gprVect))
specVect<-numeric(length(gprVect))
for(i in 1:length(gprVect)){
  sensVect[i]<-gpr95ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-gpr95ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  gpr=gprVect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-gprVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nGPR threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for GPR with diagnosis threshold te > 9.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on GPR (from 0 to 13.65 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

gpr95ResultsList[[paste(sep="","gpr",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) GPR model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 9: Best (according to Youden’s J) GPR model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 1.5598 (0.1637, 3.8903)
hazardous alcohol drinking (OR for specificity) 0.5464 (0.3089, 0.8106)
suspected liver disease screening (OR for sensitivity) 9.6354 (2.1353,20.9855)
suspected liver disease screening (OR for specificity) 0.3326 (0.2094, 0.4632)
hospital-based study (OR for sensitivity) 0.9495 (0.0844, 2.3562)
hospital-based study (OR for specificity) 0.4307 (0.3103, 0.5705)
random effects variance (logit sensitivity) 0.8032 (0.0865, 2.0909)
random effects variance (logit specificity) 0.5040 (0.1185, 1.0852)
overall sensitivity 0.8295 (0.6568, 0.9739)
overall specificity 0.7142 (0.6402, 0.7813) 
bestModels[["gpr_te95"]]<-gpr95ResultsList[[paste(sep="","gpr",bestThr)]]
bestThrs[["gpr_te95"]]<-bestThr
aucs[["gpr_te95"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="gpr",teThr=9.5)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotGprTe95",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for GPR with te >= 9.5kPa.

Figure 5: Forest plots for GPR with te >= 9.5kPa.

F2 - te > 7.9kPa

gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40

gpr79ResultsList<-list()
gpr79ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(gprVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  gpr<-gprVect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_79==1,1,2)
  datJags$y<-ifelse(dat$gpr>=gpr,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(gprVect)){
  gpr<-gprVect[i]
  
  gpr79ModList[[paste(sep="","gpr_mod_",gpr)]]<-parList[[i]]$model
  gpr79ModList[[paste(sep="","gpr_samples_",gpr)]]<-parList[[i]]$pars
  gpr79ResultsList[[paste(sep="","gpr",gpr)]]<-parList[[i]]$results
}

rm(parList)
save(gpr79ResultsList,gpr79ModList,file=paste(sep="","gprResultsList_te79_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
gprVect<-unique(round(digits=2,quantile(dat$gpr,probs=seq(0,1,length=54),na.rm=T)))
gprVect<-sort(c(0,gprVect,1,2,3,4,5,6,8,10)) # length 40

load("gprResultsList_te79_20210805.RData")
sensVect<-numeric(length(gprVect))
specVect<-numeric(length(gprVect))
for(i in 1:length(gprVect)){
  sensVect[i]<-gpr79ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-gpr79ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  gpr=gprVect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-gprVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nGPR threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for GPR with diagnosis threshold te > 7.9kPa (F2).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on GPR (from 0 to 13.65 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

gpr79ResultsList[[paste(sep="","gpr",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) GPR model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 10: Best (according to Youden’s J) GPR model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 0.7777 (0.1482, 1.6472)
hazardous alcohol drinking (OR for specificity) 0.5182 (0.2656, 0.7933)
suspected liver disease screening (OR for sensitivity) 7.0267 (2.6486,12.4961)
suspected liver disease screening (OR for specificity) 0.3135 (0.1959, 0.4376)
hospital-based study (OR for sensitivity) 1.7557 (0.3101, 3.6954)
hospital-based study (OR for specificity) 0.4872 (0.3401, 0.6424)
random effects variance (logit sensitivity) 0.8070 (0.1077, 2.0433)
random effects variance (logit specificity) 0.5383 (0.1208, 1.1604)
overall sensitivity 0.8364 (0.7037, 0.9576)
overall specificity 0.6119 (0.5305, 0.6913) 
bestModels[["gpr_te79"]]<-gpr79ResultsList[[paste(sep="","gpr",bestThr)]]
bestThrs[["gpr_te79"]]<-bestThr
aucs[["gpr_te79"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="gpr",teThr=7.9)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotGprTe79",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for GPR with te >= 7.9kPa.

Figure 6: Forest plots for GPR with te >= 7.9kPa.

Optimum thresholds for ALT

F4 - te > 11.5kPa

altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40

alt115ResultsList<-list()
alt115ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(altVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  alt<-altVect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_115==1,1,2)
  datJags$y<-ifelse(dat$alt>=alt,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(altVect)){
  alt<-altVect[i]
  
  alt115ModList[[paste(sep="","alt_mod_",alt)]]<-parList[[i]]$model
  alt115ModList[[paste(sep="","alt_samples_",alt)]]<-parList[[i]]$pars
  alt115ResultsList[[paste(sep="","alt",alt)]]<-parList[[i]]$results
}

rm(parList)
save(alt115ResultsList,alt115ModList,file=paste(sep="","altResultsList_te115_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40

load("altResultsList_te115_20210805.RData")
sensVect<-numeric(length(altVect))
specVect<-numeric(length(altVect))
for(i in 1:length(altVect)){
  sensVect[i]<-alt115ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-alt115ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  alt=altVect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-altVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nALT threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for ALT with diagnosis threshold te > 11.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on ALT (from 1 to 198 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

alt115ResultsList[[paste(sep="","alt",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) ALT model  for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 11: Best (according to Youden’s J) ALT model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 0.9139 (0.1596,1.9809)
hazardous alcohol drinking (OR for specificity) 0.7295 (0.4582,1.0217)
suspected liver disease screening (OR for sensitivity) 1.0797 (0.2004,2.2190)
suspected liver disease screening (OR for specificity) 0.2075 (0.1458,0.2730)
hospital-based study (OR for sensitivity) 1.9812 (0.1984,4.9471)
hospital-based study (OR for specificity) 0.6845 (0.5068,0.8671)
random effects variance (logit sensitivity) 0.4747 (0.0666,1.1553)
random effects variance (logit specificity) 0.8942 (0.2354,1.8993)
overall sensitivity 0.7080 (0.4767,0.9099)
overall specificity 0.7311 (0.6768,0.7824) 
bestModels[["alt_te115"]]<-alt115ResultsList[[paste(sep="","alt",bestThr)]]
bestThrs[["alt_te115"]]<-bestThr
aucs[["alt_te115"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="alt",teThr=11.5)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotAltTe115",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for ALT with te >= 11.5kPa.

Figure 7: Forest plots for ALT with te >= 11.5kPa.

F4 - te > 9.5kPa

altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40

alt95ResultsList<-list()
alt95ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(altVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  alt<-altVect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_95==1,1,2)
  datJags$y<-ifelse(dat$alt>=alt,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(altVect)){
  alt<-altVect[i]
  
  alt95ModList[[paste(sep="","alt_mod_",alt)]]<-parList[[i]]$model
  alt95ModList[[paste(sep="","alt_samples_",alt)]]<-parList[[i]]$pars
  alt95ResultsList[[paste(sep="","alt",alt)]]<-parList[[i]]$results
}

rm(parList)
save(alt95ResultsList,alt95ModList,file=paste(sep="","altResultsList_te95_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40

load("altResultsList_te95_20210805.RData")
sensVect<-numeric(length(altVect))
specVect<-numeric(length(altVect))
for(i in 1:length(altVect)){
  sensVect[i]<-alt95ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-alt95ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  alt=altVect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-altVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nALT threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for ALT with diagnosis threshold te > 9.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on ALT (from 1 to 198 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

alt95ResultsList[[paste(sep="","alt",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) ALT model  for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 12: Best (according to Youden’s J) ALT model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 1.1684 (0.3032,2.3067)
hazardous alcohol drinking (OR for specificity) 1.0205 (0.5890,1.4774)
suspected liver disease screening (OR for sensitivity) 1.7509 (0.7774,2.8940)
suspected liver disease screening (OR for specificity) 0.1432 (0.0992,0.1885)
hospital-based study (OR for sensitivity) 1.0544 (0.3193,1.9851)
hospital-based study (OR for specificity) 0.7324 (0.5285,0.9637)
random effects variance (logit sensitivity) 0.3464 (0.0635,0.8222)
random effects variance (logit specificity) 0.9517 (0.2286,2.0164)
overall sensitivity 0.5919 (0.4103,0.7667)
overall specificity 0.8213 (0.7763,0.8652) 
bestModels[["alt_te95"]]<-alt95ResultsList[[paste(sep="","alt",bestThr)]]
bestThrs[["alt_te95"]]<-bestThr
aucs[["alt_te95"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="alt",teThr=9.5)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotAltTe95",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for ALT with te >= 9.5kPa.

Figure 8: Forest plots for ALT with te >= 9.5kPa.

F2 - te > 7.9kPa

altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40

alt79ResultsList<-list()
alt79ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(altVect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  alt<-altVect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_79==1,1,2)
  datJags$y<-ifelse(dat$alt>=alt,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(altVect)){
  alt<-altVect[i]
  
  alt79ModList[[paste(sep="","alt_mod_",alt)]]<-parList[[i]]$model
  alt79ModList[[paste(sep="","alt_samples_",alt)]]<-parList[[i]]$pars
  alt79ResultsList[[paste(sep="","alt",alt)]]<-parList[[i]]$results
}

rm(parList)
save(alt79ResultsList,alt79ModList,file=paste(sep="","altResultsList_te79_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
altVect<-unique(round(digits=0,quantile(dat$alt,probs=seq(0,1,length=52))))
altVect<-c(0,altVect) # length 40

load("altResultsList_te79_20210805.RData")
sensVect<-numeric(length(altVect))
specVect<-numeric(length(altVect))
for(i in 1:length(altVect)){
  sensVect[i]<-alt79ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-alt79ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  alt=altVect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-altVect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nALT threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for ALT with diagnosis threshold te > 7.9kPa (F2).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on ALT (from 1 to 198 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

alt79ResultsList[[paste(sep="","alt",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) ALT model  for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 13: Best (according to Youden’s J) ALT model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 0.8833 (0.3293,1.5962)
hazardous alcohol drinking (OR for specificity) 0.8782 (0.4385,1.4023)
suspected liver disease screening (OR for sensitivity) 5.3853 (2.6790,8.7592)
suspected liver disease screening (OR for specificity) 0.0557 (0.0352,0.0790)
hospital-based study (OR for sensitivity) 0.9747 (0.3977,1.6336)
hospital-based study (OR for specificity) 0.6983 (0.4467,0.9820)
random effects variance (logit sensitivity) 0.8167 (0.1439,1.8155)
random effects variance (logit specificity) 1.2736 (0.2933,2.8197)
overall sensitivity 0.3709 (0.2283,0.5160)
overall specificity 0.9292 (0.9019,0.9543) 
bestModels[["alt_te79"]]<-alt79ResultsList[[paste(sep="","alt",bestThr)]]
bestThrs[["alt_te79"]]<-bestThr
aucs[["alt_te79"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="alt",teThr=7.9)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotAltTe79",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for ALT with te >= 7.9kPa.

Figure 9: Forest plots for ALT with te >= 7.9kPa.

Optimum thresholds for FIB4

F4 - te > 11.5kPa

fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43

fib4115ResultsList<-list()
fib4115ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(fib4Vect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  fib4<-fib4Vect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_115==1,1,2)
  datJags$y<-ifelse(dat$fib4>=fib4,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(fib4Vect)){
  fib4<-fib4Vect[i]
  
  fib4115ModList[[paste(sep="","fib4_mod_",fib4)]]<-parList[[i]]$model
  fib4115ModList[[paste(sep="","fib4_samples_",fib4)]]<-parList[[i]]$pars
  fib4115ResultsList[[paste(sep="","fib4",fib4)]]<-parList[[i]]$results
}

rm(parList)
save(fib4115ResultsList,fib4115ModList,file=paste(sep="","fib4ResultsList_te115_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43

load("fib4ResultsList_te115_20210805.RData")
sensVect<-numeric(length(fib4Vect))
specVect<-numeric(length(fib4Vect))
for(i in 1:length(fib4Vect)){
  sensVect[i]<-fib4115ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-fib4115ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  fib4=fib4Vect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-fib4Vect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nFIB4 threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for FIB4 with diagnosis threshold te > 11.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on FIB4 (from 0 to 38.35 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

fib4115ResultsList[[paste(sep="","fib4",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) FIB4 model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 14: Best (according to Youden’s J) FIB4 model for F4 (te > 11.5kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 1.0804 (0.1749, 2.3881)
hazardous alcohol drinking (OR for specificity) 0.6838 (0.3642, 1.0495)
suspected liver disease screening (OR for sensitivity) 5.4163 (0.9884,12.3368)
suspected liver disease screening (OR for specificity) 0.2333 (0.1094, 0.3800)
hospital-based study (OR for sensitivity) 0.8995 (0.1252, 2.0799)
hospital-based study (OR for specificity) 2.2552 (1.4912, 3.0942)
random effects variance (logit sensitivity) 1.0580 (0.1072, 2.7275)
random effects variance (logit specificity) 0.6531 (0.1213, 1.4914)
overall sensitivity 0.5790 (0.3201, 0.8425)
overall specificity 0.8794 (0.8378, 0.9191) 
bestModels[["fib4_te115"]]<-fib4115ResultsList[[paste(sep="","fib4",bestThr)]]
bestThrs[["fib4_te115"]]<-bestThr
aucs[["fib4_te115"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="fib4",teThr=11.5)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotFib4Te115",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for FIB4 with te >= 11.5kPa.

Figure 10: Forest plots for FIB4 with te >= 11.5kPa.

F4 - te > 9.5kPa

fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43

fib495ResultsList<-list()
fib495ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(fib4Vect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  fib4<-fib4Vect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_95==1,1,2)
  datJags$y<-ifelse(dat$fib4>=fib4,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(fib4Vect)){
  fib4<-fib4Vect[i]
  
  fib495ModList[[paste(sep="","fib4_mod_",fib4)]]<-parList[[i]]$model
  fib495ModList[[paste(sep="","fib4_samples_",fib4)]]<-parList[[i]]$pars
  fib495ResultsList[[paste(sep="","fib4",fib4)]]<-parList[[i]]$results
}

rm(parList)
save(fib495ResultsList,fib495ModList,file=paste(sep="","fib4ResultsList_te95_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43

load("fib4ResultsList_te95_20210805.RData")
sensVect<-numeric(length(fib4Vect))
specVect<-numeric(length(fib4Vect))
for(i in 1:length(fib4Vect)){
  sensVect[i]<-fib495ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-fib495ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  fib4=fib4Vect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-fib4Vect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nFIB4 threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for FIB4 with diagnosis threshold te > 9.5kPa (F4).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on FIB4 (from 0 to 38.35 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

fib495ResultsList[[paste(sep="","fib4",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) FIB4 model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 15: Best (according to Youden’s J) FIB4 model for F4 (te > 9.5kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 0.7012 (0.1610,1.4410)
hazardous alcohol drinking (OR for specificity) 0.7515 (0.4197,1.1372)
suspected liver disease screening (OR for sensitivity) 4.7230 (1.7261,8.3516)
suspected liver disease screening (OR for specificity) 0.3127 (0.1681,0.4883)
hospital-based study (OR for sensitivity) 0.7404 (0.1755,1.4538)
hospital-based study (OR for specificity) 2.5483 (1.7757,3.3484)
random effects variance (logit sensitivity) 0.5324 (0.0958,1.2093)
random effects variance (logit specificity) 0.3937 (0.0931,0.8414)
overall sensitivity 0.6035 (0.4076,0.7976)
overall specificity 0.7950 (0.7442,0.8447) 
bestModels[["fib4_te95"]]<-fib495ResultsList[[paste(sep="","fib4",bestThr)]]
bestThrs[["fib4_te95"]]<-bestThr
aucs[["fib4_te95"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="fib4",teThr=9.5)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotFib4Te95",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for FIB4 with te >= 9.5kPa.

Figure 11: Forest plots for FIB4 with te >= 9.5kPa.

F2 - te > 7.9kPa

fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43

fib479ResultsList<-list()
fib479ModList<-list()

cl <- parallel::makeCluster(6)
doParallel::registerDoParallel(cl)

parList<-foreach(i=1:length(fib4Vect),.packages=c("rjags","HDInterval","dplyr")) %dopar% {
  fib4<-fib4Vect[i]
  
  datJags<-list()
  
  datJags$N<-nrow(dat)
  datJags$m<-max(as.integer(factor(dat$Sitecountry)))
  datJags$study<-as.integer(factor(dat$Sitecountry))
  datJags$state<-ifelse(dat$te_79==1,1,2)
  datJags$y<-ifelse(dat$fib4>=fib4,1,0)
  datJags$studyTypeHospital<-dat$studyTypeHospital
  datJags$testReasonL<-dat$testReasonL
  datJags$alc<-dat$alc
  
  randSeed<-20210615
  set.seed(randSeed)
  
  jagsModel4 <- rjags::jags.model('jagsModel4_IndivAndStudyCovariatesNoBetweenStudyEffects_withAlc_simplerTestReason.jags', data=datJags, n.chains = 4, n.adapt = 2500)
  parsModel4 <- rjags::coda.samples(model=jagsModel4,variable.names=c('sens','spec','orGammaL','orGammaA','orLambda','reSens','reSpec'),n.iter = 6000)
  
  resultsModel4<-try(summary(parsModel4)$statistics,silent=TRUE)
  resultsModel4<-cbind(resultsModel4,t(hdi(parsModel4)))
  
  resultsModel4 <- resultsModel4 %>%
    as.data.frame() %>%
    mutate(
      Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
      upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
    mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
    dplyr::select(c(Mean,CrI))
  rownames(resultsModel4)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

  list(model=jagsModel4,pars=parsModel4,results=resultsModel4)
}

parallel::stopCluster(cl)

for(i in 1:length(fib4Vect)){
  fib4<-fib4Vect[i]
  
  fib479ModList[[paste(sep="","fib4_mod_",fib4)]]<-parList[[i]]$model
  fib479ModList[[paste(sep="","fib4_samples_",fib4)]]<-parList[[i]]$pars
  fib479ResultsList[[paste(sep="","fib4",fib4)]]<-parList[[i]]$results
}

rm(parList)
save(fib479ResultsList,fib479ModList,file=paste(sep="","fib4ResultsList_te79_",gsub(pattern="-",replacement="",Sys.Date()),".RData"))
fib4Vect<-unique(round(digits=2,quantile(dat$fib4,probs=seq(0,1,length=40))))
fib4Vect<-sort(c(0,fib4Vect,6,10)) # length 43

load("fib4ResultsList_te79_20210805.RData")
sensVect<-numeric(length(fib4Vect))
specVect<-numeric(length(fib4Vect))
for(i in 1:length(fib4Vect)){
  sensVect[i]<-fib479ResultsList[[i]]["overall sensitivity","Mean"]
  specVect[i]<-fib479ResultsList[[i]]["overall specificity","Mean"]
}
youdenJ<-sensVect+specVect-1

aucSensVect<-c(1,sensVect,0)
aucOnemSpecVect<-c(1,1-specVect,0)
idxOrder<-order(decreasing = FALSE,aucOnemSpecVect)
aucOnemSpecVect<-aucOnemSpecVect[idxOrder]
aucSensVect<-aucSensVect[idxOrder]

aucX<-aucOnemSpecVect[-1]-aucOnemSpecVect[-length(aucOnemSpecVect)]
aucY<-(aucSensVect[-length(aucSensVect)]+aucSensVect[-1])/2
auc<-sum(aucX*aucY)

dfROC<-data.frame(
  fib4=fib4Vect,
  sensitivity=sensVect,
  specificity=specVect,
  YoudenJ=youdenJ
)

idxBestYouden<-which(youdenJ==max(youdenJ))
bestThr<-fib4Vect[idxBestYouden]
bestSens<-sensVect[idxBestYouden]
bestSpec<-specVect[idxBestYouden]
dfROC %>%
  ggplot(mapping=aes(x=1-specificity,y=sensitivity)) +
  geom_smooth(col="orange",fill="orange",method="scam",formula = y ~ s(x, k = 5, bs = "mpi")) +
  geom_point(alpha=1,col="darkgrey") +
  geom_line(alpha=0.65,lwd=0.3,col="darkgrey") +
  geom_vline(xintercept=1-bestSpec,lty=2,col="grey") +
  geom_hline(yintercept=bestSens,lty=2,col="grey") +
  geom_text(x=1-bestSpec,y=bestSens,label=paste(sep="","Optimum Youden J = ",format(nsmall=2,round(digits=2,max(youdenJ))),"\nFIB4 threshold = ",bestThr)) +
  labs(title=paste(sep="","ROC curve for FIB4 with diagnosis threshold te > 7.9kPa (F2).\nAUC = ",format(nsmall=2,round(digits=2,auc))),caption="Bayesian bivariate random-effects models fitted for different threshold on FIB4 (from 0 to 38.35 according to 40 equally spaced quantiles).\nGiven the Bayesian nature of the data, there is some variability due to MCMC sampling.\nThe ROC curve (orange) is a shape-constrained generalized additive model fitted to the raw estimates (grey).\nAUC is computed from the raw estimates.") +
  theme_light()

fib479ResultsList[[paste(sep="","fib4",bestThr)]] %>%
  kable(caption="Best (according to Youden's J) FIB4 model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 16: Best (according to Youden’s J) FIB4 model for F2 (te > 7.9kPa) with site/study specific sensitivities and specificities (random effects).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 0.4472 (0.1478,0.8085)
hazardous alcohol drinking (OR for specificity) 0.9615 (0.5437,1.4109)
suspected liver disease screening (OR for sensitivity) 5.3517 (2.7271,8.4217)
suspected liver disease screening (OR for specificity) 0.4250 (0.2670,0.6050)
hospital-based study (OR for sensitivity) 0.4106 (0.1548,0.7384)
hospital-based study (OR for specificity) 1.6224 (1.1971,2.0830)
random effects variance (logit sensitivity) 0.3803 (0.0726,0.8606)
random effects variance (logit specificity) 0.4677 (0.1203,0.9764)
overall sensitivity 0.7818 (0.6616,0.8943)
overall specificity 0.5663 (0.5007,0.6301) 
bestModels[["fib4_te79"]]<-fib479ResultsList[[paste(sep="","fib4",bestThr)]]
bestThrs[["fib4_te79"]]<-bestThr
aucs[["fib4_te79"]]<-auc
doForest(dat=dat,bestThr=bestThr,marker="fib4",teThr=7.9)
## quartz_off_screen 
##                 2
include_graphics(paste(sep="","forestPlotFib4Te79",gsub(pattern="-",replacement="",Sys.Date()),".png"))
Forest plots for FIB4 with te >= 7.9kPa.

Figure 12: Forest plots for FIB4 with te >= 7.9kPa.

Appendix - models for fixed thresholds on APRI and GPR

tmp<-mod_te115_apri20 %>%
  as.data.frame() %>%
  mutate(
    Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
    upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
  mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
  dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

tmp %>%
  kable(caption="Bivariate random effects model for APRI$\\geq2.0$ and te>11.5 (F4).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 17: Bivariate random effects model for APRI\(\geq2.0\) and te>11.5 (F4).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 1.7432 (0.1091, 4.3214)
hazardous alcohol drinking (OR for specificity) Inf (0.1888, Inf)
suspected liver disease screening (OR for sensitivity) 5.2212 (0.5415,13.3010)
suspected liver disease screening (OR for specificity) 0.2905 (0.0032, 0.8721)
hospital-based study (OR for sensitivity) 0.5133 (0.0213, 1.4310)
hospital-based study (OR for specificity) 7.9510 (0.3947,24.2266)
random effects variance (logit sensitivity) 1.2211 (0.1405, 3.1288)
random effects variance (logit specificity) 1.3605 (0.0544, 4.5868)
overall sensitivity 0.2876 (0.0570, 0.5483)
overall specificity 0.9918 (0.9820, 0.9995) 
tmp<-mod_te95_apri20 %>%
  as.data.frame() %>%
  mutate(
    Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
    upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
  mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
  dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

tmp %>%
  kable(caption="Bivariate random effects model for APRI$\\geq2.0$ and te>9.5 (F4).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 18: Bivariate random effects model for APRI\(\geq2.0\) and te>9.5 (F4).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 1.7008 (0.1402, 4.0561)
hazardous alcohol drinking (OR for specificity) Inf (0.1927, Inf)
suspected liver disease screening (OR for sensitivity) 9.0142 (0.9858,22.4241)
suspected liver disease screening (OR for specificity) 0.6857 (0.0016, 2.1655)
hospital-based study (OR for sensitivity) 0.5460 (0.0182, 1.3568)
hospital-based study (OR for specificity) 8.4439 (0.2685,26.0366)
random effects variance (logit sensitivity) 1.2412 (0.1350, 3.1064)
random effects variance (logit specificity) 0.7651 (0.0649, 2.2119)
overall sensitivity 0.1476 (0.0270, 0.2852)
overall specificity 0.9916 (0.9818, 0.9990) 
tmp<-mod_te79_apri15 %>%
  as.data.frame() %>%
  mutate(
    Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
    upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
  mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
  dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

tmp %>%
  kable(caption="Bivariate random effects model for APRI$\\geq1.5$ and te>7.9 (F2).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 19: Bivariate random effects model for APRI\(\geq1.5\) and te>7.9 (F2).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 0.9711 (0.1073, 2.2673)
hazardous alcohol drinking (OR for specificity) 10.8563 (0.0749,21.1288)
suspected liver disease screening (OR for sensitivity) 18.3817 (2.7344,44.9702)
suspected liver disease screening (OR for specificity) 0.1231 (0.0050, 0.3097)
hospital-based study (OR for sensitivity) 0.8613 (0.1236, 1.8966)
hospital-based study (OR for specificity) 5.9155 (0.3831,16.7334)
random effects variance (logit sensitivity) 1.3829 (0.1927, 3.2764)
random effects variance (logit specificity) 0.8537 (0.0486, 2.5903)
overall sensitivity 0.0732 (0.0187, 0.1393)
overall specificity 0.9916 (0.9817, 0.9993) 
tmp<-mod_te115_gpr56 %>%
  as.data.frame() %>%
  mutate(
    Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
    upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
  mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
  dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

tmp %>%
  kable(caption="Bivariate random effects model for GPR$\\geq0.56$ and te>11.5 (F4).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 20: Bivariate random effects model for GPR\(\geq0.56\) and te>11.5 (F4).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 2.0999 (0.1663,5.1447)
hazardous alcohol drinking (OR for specificity) 0.5082 (0.1158,1.1127)
suspected liver disease screening (OR for sensitivity) 3.1302 (0.2450,7.9593)
suspected liver disease screening (OR for specificity) 0.3137 (0.0496,0.7137)
hospital-based study (OR for sensitivity) 2.7554 (0.1768,7.5154)
hospital-based study (OR for specificity) 1.0458 (0.3392,1.9581)
random effects variance (logit sensitivity) 1.8676 (0.1147,5.7481)
random effects variance (logit specificity) 1.0584 (0.0847,2.9765)
overall sensitivity 0.2268 (0.0333,0.4731)
overall specificity 0.9845 (0.9711,0.9956) 
tmp<-mod_te95_gpr56 %>%
  as.data.frame() %>%
  mutate(
    Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
    upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
  mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
  dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

tmp %>%
  kable(caption="Bivariate random effects model for GPR$\\geq0.56$ and te>9.5 (F4).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 21: Bivariate random effects model for GPR\(\geq0.56\) and te>9.5 (F4).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 1.6621 (0.1627, 3.8664)
hazardous alcohol drinking (OR for specificity) 0.3858 (0.0883, 0.8570)
suspected liver disease screening (OR for sensitivity) 4.8646 (0.8795,10.8084)
suspected liver disease screening (OR for specificity) 0.3742 (0.0331, 0.9331)
hospital-based study (OR for sensitivity) 1.3581 (0.2415, 3.0542)
hospital-based study (OR for specificity) 0.8700 (0.2400, 1.6434)
random effects variance (logit sensitivity) 0.8703 (0.0920, 2.3205)
random effects variance (logit specificity) 0.8122 (0.0805, 2.2122)
overall sensitivity 0.2043 (0.0534, 0.3795)
overall specificity 0.9888 (0.9788, 0.9969) 
tmp<-mod_te79_gpr32 %>%
  as.data.frame() %>%
  mutate(
    Mean=case_when(Mean>1e100~Inf,TRUE~Mean),
    upper=case_when(upper>1e100~Inf,TRUE~upper)) %>%
  mutate(CrI=paste(sep="","(",format(nsmall=4,round(digits=4,lower)),",",format(nsmall=4,round(digits=4,upper)),")")) %>%
  dplyr::select(c(Mean,CrI))
rownames(tmp)<-c("hazardous alcohol drinking (OR for sensitivity)","hazardous alcohol drinking (OR for specificity)","suspected liver disease screening (OR for sensitivity)","suspected liver disease screening (OR for specificity)","hospital-based study (OR for sensitivity)","hospital-based study (OR for specificity)","random effects variance (logit sensitivity)","random effects variance (logit specificity)","overall sensitivity","overall specificity")

tmp %>%
  kable(caption="Bivariate random effects model for GPR$\\geq0.32$ and te>7.9 (F2).",digits=4,col.names=c("Posterior mean","95% HDI credible interval"),align=c("r","r")) %>%
  kable_styling(full_width = FALSE)
Table 22: Bivariate random effects model for GPR\(\geq0.32\) and te>7.9 (F2).
Posterior mean 95% HDI credible interval
hazardous alcohol drinking (OR for sensitivity) 1.6963 (0.4509, 3.3217)
hazardous alcohol drinking (OR for specificity) 0.2790 (0.1215, 0.4660)
suspected liver disease screening (OR for sensitivity) 6.1484 (2.3796,10.9575)
suspected liver disease screening (OR for specificity) 0.3898 (0.1239, 0.7457)
hospital-based study (OR for sensitivity) 2.4301 (0.8859, 4.4199)
hospital-based study (OR for specificity) 0.8137 (0.4056, 1.2840)
random effects variance (logit sensitivity) 0.6234 (0.0888, 1.5131)
random effects variance (logit specificity) 0.8072 (0.0976, 2.0004)
overall sensitivity 0.1888 (0.0812, 0.3073)
overall specificity 0.9676 (0.9483, 0.9850)

Overall summaries

AUCs

teVect<-c("te115","te95","te79")
markerVect<-c("apri","gpr","alt","fib4")

aucDf<-data.frame(te115=rep(NA,4),te95=NA,te79=NA)
rownames(aucDf)<-toupper(markerVect)

for(te in teVect){
  for(m in markerVect){
    aucDf[toupper(m),te]<-aucs[[paste(sep="_",m,te)]]
  }
}

aucDf %>%
  kable(caption="Estimated AUCs for the different biomarkers and thresholds on transient elastrography.",digits = 2,col.names=c("te>11.5 (F4)","te>9.5 (F4)","te>7.9 (F2)")) %>%
  kable_styling(full_width = FALSE)
Table 23: Estimated AUCs for the different biomarkers and thresholds on transient elastrography.
te>11.5 (F4) te>9.5 (F4) te>7.9 (F2)
APRI 0.80 0.76 0.75
GPR 0.86 0.84 0.77
ALT 0.76 0.74 0.66
FIB4 0.80 0.74 0.72

Optimal thresholds

teVect<-c("te115","te95","te79")
markerVect<-c("apri","gpr","alt","fib4")

thrDf<-data.frame(te115=rep(NA,4),te95=NA,te79=NA)
rownames(thrDf)<-toupper(markerVect)

for(te in teVect){
  for(m in markerVect){
    thrDf[toupper(m),te]<-bestThrs[[paste(sep="_",m,te)]]
  }
}

thrDf["ALT",]<-as.character(thrDf["ALT",])

thrDf %>%
  kable(caption="Estimated optimal thresholds for the different biomarkers and thresholds on transient elastrography.",col.names=c("te>11.5 (F4)","te>9.5 (F4)","te>7.9 (F2)"),align="c") %>%
  kable_styling(full_width = FALSE)
Table 24: Estimated optimal thresholds for the different biomarkers and thresholds on transient elastrography.
te>11.5 (F4) te>9.5 (F4) te>7.9 (F2)
APRI 0.53 0.44 0.44
GPR 0.13 0.13 0.11
ALT 27 31 38
FIB4 1.65 1.44 0.98

Sensitivities (for the optimal threshold models)

teVect<-c("te115","te95","te79")
markerVect<-c("apri","gpr","alt","fib4")

sensDf<-data.frame(te115=rep(NA,4),te95=NA,te79=NA)
rownames(sensDf)<-toupper(markerVect)

for(te in teVect){
  for(m in markerVect){
    sensDf[toupper(m),te]<-paste(sep=" ",format(nsmall=4,round(digits=4,bestModels[[paste(sep="_",m,te)]]["overall sensitivity","Mean"])),bestModels[[paste(sep="_",m,te)]]["overall sensitivity","CrI"])
  }
}

sensDf %>%
  kable(caption="Estimated sensitivities (95% confidence intervals) for the different biomarkers and thresholds on transient elastrography.",digits = 2,col.names=c("te>11.5 (F4)","te>9.5 (F4)","te>7.9 (F2)")) %>%
  kable_styling(full_width = FALSE)
Table 25: Estimated sensitivities (95% confidence intervals) for the different biomarkers and thresholds on transient elastrography.
te>11.5 (F4) te>9.5 (F4) te>7.9 (F2)
APRI 0.6222 (0.3468,0.8748) 0.7049 (0.5192, 0.8895) 0.6674 (0.5289, 0.8061)
GPR 0.8622 (0.6569, 0.9994) 0.8295 (0.6568, 0.9739) 0.8364 (0.7037, 0.9576)
ALT 0.7080 (0.4767,0.9099) 0.5919 (0.4103,0.7667) 0.3709 (0.2283,0.5160)
FIB4 0.5790 (0.3201, 0.8425) 0.6035 (0.4076,0.7976) 0.7818 (0.6616,0.8943)

Specificities (for the optimal threshold models)

teVect<-c("te115","te95","te79")
markerVect<-c("apri","gpr","alt","fib4")

specDf<-data.frame(te115=rep(NA,4),te95=NA,te79=NA)
rownames(specDf)<-toupper(markerVect)

for(te in teVect){
  for(m in markerVect){
    specDf[toupper(m),te]<-paste(sep=" ",format(nsmall=4,round(digits=4,bestModels[[paste(sep="_",m,te)]]["overall specificity","Mean"])),bestModels[[paste(sep="_",m,te)]]["overall specificity","CrI"])
  }
}

specDf %>%
  kable(caption="Estimated specificities (95% confidence intervals) for the different biomarkers and thresholds on transient elastrography.",digits = 2,col.names=c("te>11.5 (F4)","te>9.5 (F4)","te>7.9 (F2)")) %>%
  kable_styling(full_width = FALSE)
Table 26: Estimated specificities (95% confidence intervals) for the different biomarkers and thresholds on transient elastrography.
te>11.5 (F4) te>9.5 (F4) te>7.9 (F2)
APRI 0.8564 (0.8162,0.8962) 0.7750 (0.7263, 0.8235) 0.8091 (0.7597, 0.8563)
GPR 0.6979 (0.6299, 0.7637) 0.7142 (0.6402, 0.7813) 0.6119 (0.5305, 0.6913)
ALT 0.7311 (0.6768,0.7824) 0.8213 (0.7763,0.8652) 0.9292 (0.9019,0.9543)
FIB4 0.8794 (0.8378, 0.9191) 0.7950 (0.7442,0.8447) 0.5663 (0.5007,0.6301)